Nested Grassmannians for Dimensionality Reduction with Applications

In the recent past, nested structures in Riemannian manifolds has been studied in the context of dimensionality reduction as an alternative to the popular principal geodesic analysis (PGA) technique, for example, the principal nested spheres. In this paper, we propose a novel framework for constructing a nested sequence of homogeneous Riemannian manifolds. Common examples of homogeneous Riemannian manifolds include the n-sphere, the Stiefel manifold, the Grassmann manifold and many others. In particular, we focus on applying the proposed framework to the Grassmann manifold, giving rise to the nested Grassmannians (NG). An important application in which Grassmann manifolds are encountered is planar shape analysis. Specifically, each planar (2D) shape can be represented as a point in the complex projective space which is a complex Grassmann manifold. Some salient features of our framework are: (i) it explicitly exploits the geometry of the homogeneous Riemannian manifolds and (ii) the nested lower-dimensional submanifolds need not be geodesic. With the proposed NG structure, we develop algorithms for the supervised and unsupervised dimensionality reduction problems respectively. The proposed algorithms are compared with PGA via simulation studies and real data experiments and are shown to achieve a higher ratio of expressed variance compared to PGA.


Introduction
Riemannian manifolds are often used to model the sample space in which features derived from the raw data encountered in many medical imaging applications live. examples include the diffusion tensors (DTs) in diffusion tensor imaging (DTI) (Basser et al., 1994), the ensemble average propagator (EAP) (Callaghan, 1993). Both DTs and EAP are used to capture the diffusional properties of water molecules in the central nervous system by non-invasively imaging the tissue via diffusion weighted magnetic resonance imaging. In DTI, diffusion at a voxel is captured by a DT, which is a 3 × 3 symmetric positive-definite matrix, whereas EAP is a probability distribution characterizing the local diffusion at a voxel, which can be parametrized as a point on the Hilbert sphere. Another example is the shape space used to represent shapes in shape analysis. There are many ways to represent a shape, and the most simple one is to use landmarks. For the landmark-based representation, the shape space is called Kendall's shape space (Kendall, 1984). Kendall's shape space is in general a stratified space (Goresky and MacPherson, 1988;Feragen et al., 2014), but for the special case of planar shapes, the shape space is the complex projective space, which is a complex Grassmann manifold. The examples mentioned above are often high-dimensional: a DTI scan usually contains half a million DTs; the shape of the Corpus Callosum (which is used in our experiments) is represented by a several hundreds of boundary points in ℝ 2 . Thus, in these cases, dimension reduction techniques, if applied appropriately, can benefit the subsequent statistical analysis.
For data on Riemannian manifolds, the most widely used dimensionality reduction method is the principal geodesic analysis (PGA) (Fletcher et al., 2004), which generalizes the principal component analysis (PCA) to manifold-valued data. In fact, there are many variants of PGA. Fletcher et al. (2004) proposed to find the geodesic submanifold of a certain dimension that maximizes the projected variance and computationally, it was achieved via a linear approximation, i.e., applying PCA on the tangent space at the intrinsic mean. This is sometimes referred to as the tangent PCA. Note that this approximation requires the data to be clustered around the intrinsic mean, otherwise the tangent space approximation to the manifold leads to inaccuracies. Later on, Sommer et al. (2010) proposed the Exact PGA (EPGA), which does not use any linear approximation. However, EPGA is computationally expensive as it requires two non-linear optimizations steps per iteration (projection to the geodesic submanifold and finding the new geodesic direction such that the loss of information is minimized). Chakraborty et al. (2016) partially solved this problem for manifolds with constant sectional curvature (spheres and hyperbolic spaces) by deriving closed form formulae for the projection. Other variants of PGA include but are not limited to sparse exact PGA (Banerjee et al., 2017), geodesic PCA (Huckemann et al., 2010), and probabilistic PGA (Zhang and Fletcher, 2013). All these methods focus on projecting data to a geodesic submanifold as in PCA where one projects data to a vector subspace. Instead, one can also project data to a submanifold that minimizes the reconstruction error without any further restrictions, e.g. being geodesic. This is the generalization of the principal curve (Hastie and Stuetzle, 1989) to Riemannian manifolds presented in Hauberg (2016). the manifold P n of n × n SPD matrices, Harandi et al. (2018) proposed a geometry-aware dimension reduction by projecting data in P n to P m for some m ≪ n. They also applied the nested P n for the supervised dimensionality reduction problem. Damon and Marron (2014) considered a nested sequence of relations which determine a nested sequence of submanifolds that are not necessarily geodesic. They showed various examples, including Euclidean space and the n-sphere, depicting how the nested relations generalized PCA and PNS. However, for an arbitrary Riemannian manifold, it is not clear how to construct a nested submanifold. Another generalization of PGA was proposed by Pennec et al. (2018), called the exponential barycentric subspace (EBS). A k-dimensional EBS is defined as the locus of weighted exponential barycenters of (k + 1) affinely independent reference points. The EBSs are naturally nested by removing or adding reference points.
Unlike PGA which can be applied to any Riemannian manifolds, the construction of the nested manifolds relies heavily on the geometry of the specific manifold, and there is no general principle for such a construction. All the examples described above (spheres and P n ) and several others such as the Grassmannian, Stiefel etc. are homogeneous Riemannian manifolds (Helgason, 1979), which intuitively means that all points on the manifold 'look' the same. In this work, we propose a general framework for constructing a nested sequence of homogeneous Riemannian manifolds, and, via some simple algebraic computation, show that the nested sphere and the nested P n can indeed be derived within this framework. We will then apply this framework to the Grassmann manifolds, called the nested Grassmann manifolds (NG). The Grassmann manifold Gr(p, V) is the manifold of all p-dimensional subspaces of the vector space V where 1 ≤ p ≤ dim V. Usually V = ℝ n or V = ℂ n . An important example is Kendall's shape space of 2D shapes. The space of all shapes determined by k landmarks in ℝ 2 is denoted by Σ 2 k , and Kendall (1984) showed that it is isomorphic to the complex projective space ℂP k − 2 ≅ Gr 1, ℂ k − 1 . In many applications, the number k of landmarks is large, and so is the dimension of Gr 1, ℂ k − 1 . The core of the proposed dimensionality reduction involves projecting data on Gr(p, V) to Gr(p, V) with dim V ≪ dim V. The main contributions of this work are as follows: (i) We propose a general framework for constructing a nested sequence of homogeneous Riemannian manifolds unifying the recently proposed nested spheres (Jung et al., 2012) and nested SPD manifolds (Harandi et al., 2018). (ii) We present novel dimensionality reduction techniques based on the concept of NG in both supervised and unsupervised settings respectively. (iii) We demonstrate via several simulation studies and real data experiments, that the proposed NG can achieve a higher ratio of expressed variance compared to PGA.
The rest of the paper is organized as follows. In Section 2, we briefly review the definition of homogeneous Riemannian manifolds and present the recipe for the construction of nested homogeneous Riemannian manifolds. In Section 3, we first review the geometry of the Grassmannian. By applying the procedure developed in Section 2, we present the nested Grassmann manifolds representation and discuss some of its properties in details. Then we describe algorithms for our unsupervised and supervised dimensionality reduction techniques, called the Principal Nested Grassmanns (PNG), in Section 4. In Section 5, we present several simulation studies and experimental results showing how the PNG technique performs compared to PGA under different settings. Finally, we draw conclusions in Section 6.

Nested Homogeneous Spaces
In this section, we introduce the structure of nested homogeneous Riemannian manifolds. A Riemannian manifold (M, τ) is homogeneous if the group of isometries G = I(M) admitted by the manifold acts transitively on M (Helgason, 1979), i.e., for x, y ∈ M, there exists g ∈ G such that g(x) = y. In this case, M can be identified with G/H where H is an isotropy subgroup of G at some point p ∈ M i.e. H = {g ∈ G : g(p) = p}. Examples of homogeneous Riemannian manifolds include but are not limited to, the n-spheres S n−1 = SO(n)/SO(n -1), the SPD manifolds P n = GL(n)/O(n), the Stiefel manifolds St(m, n) = SO(n)/SO(n − m), and the Grassmann manifolds Gr(p, n) = SO(n)/S(O(p) × O(n − p)).
In this paper, we focus on the case where G is either a real or a complex matrix Lie group, i.e. G is a subgroup of GL(n, ℝ) or GL(n, ℂ). The main idea behind the construction of nested homogeneous spaces is simple: augmenting the matrix in G in an appropriate way. With an embedding of the isometry group G, the embedding of the homogeneous space G/H follows naturally from the quotient structure.
Let G and G be two connected Lie groups such that dim G < dim G and ι : G G be an embedding. For a closed connected subgroup H of G, let H = ι (H). Since ι is an embedding, H is also a closed subgroup of G. Now the canonical embedding of G/H in G/H is defined by ι(gH) = ι (g)H for g ∈ G. It is easy to see that ι is well-defined. Let g 1 , g 2 ∈ G be such that g 1 = g 2 h for some h ∈ H. Then ι g 1 H = ι g 1 H = ı g 2 ℎ H = ι g 2 ι (ℎ)H = ι g 2 H = ι g 2 H . Now for the homogeneous Riemannian manifolds (M = G/H, τ 1 ) and M = G/H, τ 2 , denote the left-G-invariant, right-H-invariant metric on G (resp. left-G-invariant, right-H-invariant metric on G) by τ 1 and τ 2 , respectively (see Cheeger and Ebin (1975, Prop. 3.16(4))).
Proposition 1 If ι : G G is isometric, then so is ι: G/H G/H. Proof Denote the Riemannian submersions by π : G → G/H and π: G G/H. Let X and Y be vector fields on G/H and X and Y be their horizontal lifts respectively, i.e. dπ(X) = X and dπ(Y ) = Y . By the definition of Riemannian submersions, dπ is isometric on the horizontal spaces, i.e. τ 1 (X, Y ) = τ 1 (dπ(X), dπ(Y )) = τ 1 (X, Y ). Since ı is isometric, we have τ 1 (X, Y ) = τ 2 (dι (X), dι (Y )). By the definition of ι, we also have ι ∘ π = π ∘ ι , which implies where the third equality follows from the isometry of dπ. ■ Proposition 1 simply says that if the isometry group is isometrically embedded, then the associated homogeneous Riemannian manifolds will also be isometrically embedded. Conversely, if we have a Riemannian submersion f : G G, it can easily be shown that the induced map f : G/H G/H would also be a Riemannian submersion where H = f(H). The construction above can be applied to a sequence of homogeneous spaces M m m = 1 ∞ , i.e. the embedding ι m : M m → M m+1 can be induced from the embedding of the isometry groups ι m : G m G m + 1 where G m = I(M m ) provided dim G i < dim G j for i < j. See Figure 1 for the structure of nested homogeneous spaces.

Nested Grassmann Manifolds
In this section, we will apply the theory of nested homogeneous space from the previous section to the Grassmann manifolds. First, we briefly review the geometry of the Grassmann manifolds in Section 3.1. With the theory in Section 2, we derive the nested Grassmann manifolds in Section 3.2, and the derivation for nested spheres and nested SPD manifolds are carried out in Section 3.3.

The Riemannian Geometry of Grassmann Manifolds
To simplify the notation, we assume V = ℝ n and write Gr(p, n) ≔ Gr p, ℝ n . All the results presented in this section can be easily extended to the case of V = ℂ n by replacing transposition with conjugate transposition and orthogonal groups with unitary groups. The Grassmann manifold Gr(p, n) is the manifold of all p-dimensional subspaces of ℝ n , and for a subspace ∈ Gr(p, n),. we write = span(X) where the columns of X form an orthonormal basis for . The space of all n × p matrices X such that X T X = I p called the Stiefel manifold, denoted by St(p, n). Special cases of Stiefel manifolds are the Lie group of all orthogonal matrices, O(n) = St(n, n), and the n-sphere, S n−1 = St(1, n). The Stiefel manifold with the induced Euclidean metric (i.e. for U, V ∈ T X St(p, n). 〈U, V〉 X = tr(U T V)) is a homogeneous Riemannian manifold, St(p, n) = SO(n)/SO(n − p). A canonical Riemannian metric on the Grassmann manifold can be inherited from the metric on St(p, n) since it is invariant to the left multiplication by elements of O(n) (Absil et al., 2004;Edelman et al., 1998). The Grassmann manifold with this metric is also homogeneous, Gr(p, With this canonical metric on the Grassmann manifolds, the geodesic can be expressed in closed form. Let = span(X) ∈ Gr(p, n) where X ∈ St(p, n) and H be an n × p matrix. Then the geodesic γ(t) with γ(0) = and γ′(0) = H given by γ ,H (t) = span(XV cos Σt + U sin Σt) where H = UΣV T is the compact singular value decomposition (Edelman et al., 1998, Theorem 2.3). The exponential map at is a map from T Gr(p, n) to Gr(p, n) defined by Exp H = γ ,H (1) = span(XV cos Σ + U sin Σ). If X T Y is invertible, the geodesic distance between = span(X) and = span(Y) is given by d g

Embedding of Gr(p, m) in Gr(p, n)
Let = span(X) ∈ Gr(p, m), X ∈ St(p, m). The map ι : Gr(p, m) → Gr(p, n), for m < n, defined by ι(X) = span X 0 (n − m) × p is an embedding and it is easy to check that this embedding is isometric (Ye and Lim, 2016, Eq. (8)). However, for the dimensionality reduction problem, the above embedding is insufficient as it is not flexible enough to encompass other possible embeddings. To design flexible embeddings, we apply the framework proposed in Section 2. We consider M m = Gr(p, m) for which the isometry groups are G m = SO(m) and H m = S(O(p) × O(m − p)).
In this paper, we consider the embedding ι m : SO(m) SO(m + 1) given by, and O p is the m × p matrix containing the first p columns of O, the induced embedding ι m : Gr(p, m) → Gr(p, m + 1) is given by, where b ∈ ℝ p , R ∈ SO(m + 1), R contains the first m columns of R (which means R ∈ St m, m + 1 , υ is the last column of R, and = span(X) ∈ Gr(p, m). It is easy to see that for R = I and b = 0, this gives the natural embedding described in Ye and Lim (2016) and at the beginning of this section.
Proposition 2 If b = 0, then ι m is an isometric embedding.
Proof With Proposition 1, it suffices to show that ι m is isometric when b = 0. Note that as ι m is independent of a and c in the definition of ι m , we can assume a = 0 and c = 1 without loss of generality. If ι m simplifies to where R ∈ SO(m + 1). The Riemannian distance on SO(n) given the induced Euclidean Therefore ι m is an isometric embedding, and so is ι m by Proposition 1. ■ With the embedding ι m , we can construct the corresponding projection π m : Gr(p, m + 1) → Gr(p, m) using the following proposition.
the projection is given by π m (X) = span R T X . This completes the proof. ■ Note that for = span(X) 2 Gr(p, m+1), ι m (π m ( )) = span(RR T X+vb T ) = span((I − vv T )X +vb T ) where v ∈ ℝ m + 1 and ∥υ∥ = 1. The nested relation can be extended inductively and we refer to this construction as the nested Grassmann structure: Thus the embedding from Gr(p, m) into Gr(p, n) can be constructed inductively by ι ≔ ι n−1 ∘… ∘ι m−1 ∘ ι m and similarly for the corresponding projection. The explicit forms of the embedding and the projection are given in the following proposition.

Proposition 4
The embedding of Gr(p, m) into Gr(p, n) for m < n is given by ι A,B ( ) = span(AX + B) where A ∈ St(m, n) and B ∈ ℝ n × p such that A T B = 0. The corresponding projection from Gr(p, n) to Gr(p, m) is given by π A = span(A T X).
Proof By the definition, ι ≔ ι n−1 ∘…∘ι m−1 ∘ι m and thus the embedding ι : Gr(p, m) → Gr(p, n) can be simplified as T is an n × p matrix. Similar to Proposition 3, the projection π A : Gr(p, n) → Gr(p, m) is then given by π A ( ) = span(A T X). This completes the proof. ■ From Proposition 2, if B = 0 then ι A is an isometric embedding. Hence, our nested Grassmann structure is more flexible than PGA as it allows one to project the data onto a non-geodesic submanifold. An illustration is shown in Figure 2. The results discussed in this subsection can be generalized to any homogeneous space in principle. For a given homogeneous space, once the embedding of the groups of isometries (e.g., Eq. (1)) is determined, the induced embedding and the corresponding projection can be derived akin to the case of Grassmann manifolds.

Connections to Other Nested Structures
The nested homogeneous spaces proposed in this work (see Figure 1) actually provides a unified framework within which, the nested spheres (Jung et al., 2012) and the nested SPD manifolds (Harandi et al., 2018) are special cases.
The n-Sphere Example: Since the n-sphere can be identified with a homogeneous space S n−1 ≅ O(n)/O(n − 1), with the embedding (1), the induced embedding of S n−1 into S n is where x ∈ S n−1 , b ∈ ℝ, and r = cos −1 b . This is precisely the nested sphere proposed

Dimensionality Reduction with Nested Grassmanns
In this section, we discuss how to apply the nested Grassmann structure to the problem of dimension reduction. In Section 4.1 and 4.2, we describe the unsupervised and supervised dimension reduction based on the nested Grassmann manifolds. In Section 4.3, we will discuss the choice of distance metrics required by the dimensionality reduction algorithm and present some technical details regarding the implementation. Analysis of principal nested Grassmann (PNG) will be introduced and discussed in Section 4.4 and Section 4.5.

Unsupervised Dimensionality Reduction
We can now apply the nested Grassmann structure to the problem of unsupervised dimensionality reduction. Suppose that we are given the points, that we seek is obtained by the minimizing the reconstruction error, i.e. 1 , …, N ∈ Gr(p, n). We would like to have lower dimensional representations in Gr(p, m) for 1 , …, N with m ≪ n. The desired projection map π A that we seek is obtained by the minimizing the reconstruction error, i.e. and it is optimized over Gr(m, n) × ℝ n × p . Note that since we represent a subspace by its orthonormal basis, when m > n/2, we can use the isomorphism Gr(m, n) ≅ Gr(n − m, n) to further reduce the computational burden. This will be particularly useful when m = n − 1 as in Section 4.4. Under this isomorphism Gr(m, n) ≅ Gr(n − m, n), the corresponding subspace of span (A) ∈ Gr(m, n) is span(A ⊥ ) ∈ Gr(n − m, n) where A ⊥ is an n × (n − m) matrix such that [A A ⊥ ] is an orthogonal matrix. Hence the loss function L u can alternatively be expressed as Yang

Supervised Dimensionality Reduction
If in addition to 1 , …, N ∈ Gr(p, n), we are given the associated labels y 1 , …, y N ∈ {1, …, k}, then we would like to use this extra information to sharpen the result of dimensionality reduction. Specifically, we expect that after reducing the dimension, points from the same class preserve their proximity while points from different classes are well separated. We use an affinity function a: Gr(p, n) × Gr(p, n) ℝ to encode the structure of the data as suggested by Harandi et al. (2018, Sec 3.1, Eq. (14)- (16)). For completeness, we repeat the definition of the affinity function here. The affinity function is defined as a( i , The set Nw( i ) consists of ν w nearest neighbors of i that have the same labels as y i , and the set N b ( i ) consists of ν b nearest neighbors of i that have different labels from y i . The nearest neighbors can be computed using the geodesic distance.
The desired projection map π A that we seek is obtained by the minimizing the following loss function

Choice of the distance function
The loss functions L u and L s depend on the choice of the distance d: Gr(p, n) × Gr(p, n) ℝ ≥ 0 . Besides the geodesic distance, there are many widely used distances on the Grassmann manifold, see, for example, Edelman et al. (1998, p. 337) Yang and Vemuri Page 10 and Ye and Lim (2016, Table 2). In this work, we use two different distance metrics: (1) the geodesic distance d g and (2) the projection distance, which is also called the chordal distance in Ye and Lim (2016) and the projection F-norm in Edelman et al.
(1998). The geodesic distance was defined in Section 3.1 and the projection distance is defined as follows. For , ∈ Gr(p, n), denote the projection matrices onto and is given by p and p respectively. Then, the distance between and is given by where θ 1 , …, θ p are the principal angles of and . If = span(X) then P = X(X T X) −1 X T . It is also easy to see the the projection distance has O(n)-symmetry. We choose the projection distance mainly for its computational efficiency as it involves only matrix multiplication which has a time complexity O(n 2 ) whereas the geodesic distance requires an SVD which has a time complexity of O(n 3 ).

Analysis of Principal Nested Grassmannians
To determine the dimension of the nested submanifold that fits the data well enough, we can choose p < m 1 < … < m k < n and estimate the projection onto these nested Grassmann manifolds. The ratio of expressed variance for each projection is the ratio of the variance of the projected data and the variance of the original data. With these ratios, we can choose the desired dimension according to the pre-specified percentage of expressed variance as one would do for choosing the number of principal components in PCA.
Alternatively, one can have a full analysis of principal nested Grassmanns (PNG) as follows.
Furthermore, we can reduce the points on Gr(1, 2), which is a 1-dimensional manifold, to a 0-dimensional manifold, which is simply a point, by computing the Fréchet mean (FM). We call this FM the nested Grassmannian mean (NGM) of 1 , …, N ∈ Gr(p, n). The NGM is unique since Gr(1, 2) ≅ ℝP 1 can be identified as the half circle in ℝ 2 and the FM is unique in this case. Note that in general, the NGM will not be the same as the FM of 1 , …, N since the embedding/projection need not be isometric. The supervised PNG (sPNG) can be obtained similarly by replacing each projection with it supervised counterpart.

Principal Scores
Whenever we apply a projection π m : Gr(p, m + 1) → Gr(p, m) to the data, we might lose some information contained in the data. More specifically, since we project data on a p(m + 1 − p)-dimensional manifold to a p(m − p)-dimensional manifold, we need to describe this p(m + 1 − p) p(m − p) = p dimensional information loss during the projection. In PCA, this is done by computing the scores of each principal component, which are the transformed coordinates of each sample in the eigenspace of the covariance matrix. We can generalize the notion of principal scores to the nested Grassmanns as follows: For each ∈ Gr(p, m + 1), denote by M π m (X) the fiber of π m ( ), i.e. M π m (X) = π m −1 π m (X) = Y ∈ Gr(p, m + 1): π m (Y) = π m (X) which is a p-dimensional submanifold of Gr(p, m + 1). An illustration of this fiber is given in Figure 3a. Let X = ι m π m (X) and let the unit tangent vector V ∈ T X M π m (X) be the geodesic direction from X to . Given a suitable basis on T X M π m (X) , V can be realized as a p-dimensional vector, and this will be the score vector of associated with the projection, π m.
By the definition of M π m (X) we have the following decomposition of the tangent space of Gr(p, m + 1) at X into the horizontal space and the vertical space induced by π m , T X Gr(p, m + 1) = T X M π m (X) ⊕ dι m π m (X) T π m (X) Gr(p, m) .
An illustration of this decomposition is given in Figure 3b. A tangent vector of M π m (X) at X is of the form Δ = A ⊥ b T where A ⊥ is any (m + 1)-dim vector such that [A A ⊥ ] is orthogonal and b ∈ ℝ p It is easy to check that π m (span(AA T X + A ⊥ b T )) = π m (span(X)) = span(A T X).
Hence a natural coordinate for the tangent vector Δ = A ⊥ b T is b ∈ ℝ p , and the geodesic direction from X to would be V= X T A ⊥. It is easy to see that ∥V∥ F = 1 since X has orthonormal columns. To reflect the distance between X and , i.e. the reconstruction error, we define d X, X V as the score vector for associated with π m . In the case of Gr(1, 2) → NGM, we use the sign distance to the NGM as the score. For complex nested Grassmanns however, the principal score associated with each projection is a p-dimensional complex vector. For the sake of visualization, we transform this p-dimensional complex vector to a 2p-dimensional real vector. The procedure for computing the PNG and the principal scores is summarized in Algorithm 1.
Remark 2 Note that this definition of principal score is not intrinsic as it depends on the choice of basis. Indeed, it is impossible to choose a p-dimensional vector for the projection π m in an intrinsic way, since the only property of a map that is independent of bases is the rank of the map. A reasonable choice of basis is made by viewing the Grassmann Gr(p, m) as a quotient manifold of St(p, m), which is a submanifold in ℝ m × p . This is how we define the principal score for nested Grassmanns.

Experiments
In this section, we will demonstrate the performance of the proposed dimensionality reduction technique, i.e. PNG and sPNG, via experiments on synthetic and real data. The implementation 1 is based on the python library pymanopt (Townsend et al., 2016) and we use the steepest descent algorithm for the optimization (with default parameters in pymanopt). The optimization was performed on a desktop with 3.6GHz Intel i7 processors and took about 30 seconds to converge.

Synthetic Data
In this subsection, we compare the performance of the projection and the geodesic distances respectively. The questions we will answer are the following.
(1) From Section 4.3, we see that using projection distance is more efficient than using the geodesic distance. But how do they perform compared to each other under varying dimension n and variance level σ 2 ? (2) Is our method of dimensionality reduction "better" than PGA? Under what conditions does our method outperform PGA? Finally, we compute X i = Exp X i σU i , where U i = U i / U i and U i ∈ T X i Gr(p, n) to include some perturbation.

Projection and Geodesic
This experiment involves comparing the performance of the NG representation in terms of the explained variance, under different levels of data variance. In this experiment, we set N = 50, n = 10, m = 3, and p = 1 and σ is ranging from 0.5 to 1. The results are averaged over 100 repetitions and are shown in Figure 4. From these results, we can see that the explained variance for the projection distance and the geodesic distance are indistinguishable but using projection distance leads to much faster convergence than when using the geodesic distance. The reason is that when two points on the Grassmann manifold are close, the geodesic distance can be well-approximated by the projection distance. When the algorithm converges, the original point i and the reconstructed point X i should be close and the geodesic distance can thus be well-approximated by the projection distance. Therefore, for all the experiments in the next section, we use the projection distance for the sake of efficiency.

Comparison of PNG and PGA-Now
we compare the proposed PNG to PGA. In order to have a fair comparison between PNG and PGA, we define the principal components of PNG as the principal components of the scores obtained as in Section 4.5. Similar to the previous experiment, we set N = 50, n = 10, m = 5, p = 2, and σ = 0.01, 0.05, 0.1, 0.5 and apply the same procedure to generate synthetic data. The results are averaged over 100 repetitions and are shown in Figure 5.
From Figure 5, we can see that our method outperforms PGA by virtue of the fact that it is able to capture a larger amount of variance contained in the data. We can see that when the variance is small, the improvement of PNG over PGA is less significant, whereas, our method is significantly better for the large data variance case (e.g. comparing σ = 0.5 and σ = 0.01). Note that when the variance in the data is small, i.e. the data are tightly clustered around the FM, and PGA captures the essence of the data well. However, the requirement in PGA on the geodesic submanifold to pass through the anchor point, namely the FM, is not meaningful for data with large variance as explained through the following simple example. Consider, a few data points spread out on the equator of a sphere. The FM in this case is likely to be the north pole of the sphere if we restrict ourselves to the upper hemisphere. Thus, the geodesic submanifold computed by PGA will pass through this FM. However, what is more meaningful is a submanifold corresponding to the equator, which is what a nested spheres representation (Jung et al., 2012) in this case yields. In similar vein, for data with large variance on a Grassmann manifold, our NG representation will yield a more meaningful representation than PGA.

Application to Planar Shape Analysis
We now apply our method to planar (2-dimensional) shape analysis. A planar shape σ can be represented as an ordered set of k > 2 points in ℝ 2 , called a k-ad or a configuration. Here we assume that these k points are not all identical. Denote the configuration by X which is a k × 2 matrix. Let H be the (k − 1) × k Helmert submatrix (Dryden and Mardia, 2016, Ch. 2.5). Then Z = HX/∥HX∥ F is called the pre-shape of X from which the information about translation and scaling is removed. The space of all pre-shapes is called the pre-shape space, denoted by S 2 k . By definition, the pre-shape space is a (2k − 3)-dimensional sphere.
The shape is obtained by removing the rotation from the pre-shape, and thus the shape space is Σ 2 k = S 2 k /O(2). It was shown by Kendall (1984) that Σ 2 k is a smooth manifold and, when equipped with the quotient metric, is isometric to the complex projective space ℂP k − 2 equipped with the Fubini-Study metric (up to a scale factor) which is a special case of the complex Grassmannians, i.e. ℂP k − 2 ≅ Gr 1, ℂ k − 1 . Hence, we can apply the proposed PNG to planar shapes. For planar shapes, we also compare with the recently proposed principal nested shape spaces (PNSS) (Dryden et al., 2019), which is an application of PNS on the pre-shape space. We will now demonstrate how the PNG performs compared to PGA and PNSS using some simple examples of planar shapes and the OASIS dataset.
Examples of Planar Shapes-We take three datasets: digit3, gorf, and gorm, from the R package shapes (Dryden, 2021). The digit3 dataset consists of 30 shapes of the digit 3, each of which is represented by 13 points in ℝ 2 ; the gorf dataset consists of 30 shapes of female gorilla skull, each of which is represented by 8 points in ℝ 2 ; the gorm dataset consists of 29 shapes of male gorilla skull, each of which is represented by 8 points in ℝ 2 . Example shapes from these three datasets are shown in Figure 6. The cumulative ratios of variance explained by the first 5 principal components 2 of PNG, PGA, and PNSS are shown in Figure 7. It can be seen from Figure 7 that the proposed PNG achieves higher explained variance than PGA and PNSS respectively in most cases.

OASIS Corpus Callosum Data Experiment-
The OASIS database (Marcus et al., 2007) is a publicly available database that contains T1-MR brain scans of subjects of age ranging from 18 to 96. In particular, it includes subjects that are clinically diagnosed with mild to moderate Alzheimer's disease. We further classify them into three groups: young (aged between 10 and 40), middle-aged (aged between 40 and 70), and old (aged above 70).
For demonstration, we randomly choose 4 brain scans within each decade, totalling 36 brain scans. From each scan, the Corpus Callosum (CC) region is segmented and 250 points are taken on the boundary of the CC region. See Figure 8 for samples of the segmented corpus callosi. In this case, the shape space is Σ 2 248 ≅ ℂP 248 ≅ Gr 1, ℂ 249 . Results of application of the three methods to this data are shown in Figure 9.
Since the data are divided into three groups (young, middle-aged, and old), we can apply the sPNG described in Section 4.2 to reduce the dimension. The purpose of this experiment is not to demonstrate state-of-the-art classification accuracy for this dataset. Instead, our goal here is to demonstrate that the proposed nested Grassmann representation in a supervised setting is much more discriminative than the competition, namely the supervised PGA.
Hence, we choose a simple classifier such as the support vector machine (SVM) Vapnik (1995) to highlight the aforementioned discriminative power of the nested Grassmann over PGA.
For comparison, the PGA can be easily extended to supervised PGA (sPGA) by first diffeomorphically mapping all the data to the tangent space anchored at the FM and then performing supervised PCA Bair et al. (2006); Barshan et al. (2011) on the tangent space. However, generalizing PNSS to the supervised case is nontrivial and is beyond the scope of this paper. Therefore, we limit our comparison to the unsupervised PNSS. In this demonstration, we apply an SVM on the scores obtained from different dimension reduction algorithms, and we choose only the first three principal scores to show that even with the 3-dimensional representation of the original shapes, we can still achieve good classification results. The results are shown in Table 1. These results are in accordance with our expectation since in sPNG, we seek a projection that minimizes the within-group variance while maximizing the between-group variance. However, as we observed earlier, the constraint of requiring the geodesic submanifold to pass through the FM is not well suited for this dataset which has a large variance across the data. This accounts for why the sPNG exhibits far superior performance compared to sPGA in accuracy.

Conclusion
In this work, we proposed a novel nested geometric structure for homogeneous spaces and used this structure to achieve dimensionality reduction for data residing in Grassmann manifolds. We also discuss how this nested geometric structure served as a natural generalization of other existing nested geometric structures in literature namely, spheres and the manifold of SPD matrices. Specifically, we showed that a lower dimensional Grassmann manifold can be embedded into a higher dimensional Grassmann manifold and via this embedding we constructed a sequence of nested Grassmann manifolds. Compared to the PGA, which is designed for general Riemannian manifolds, the proposed method can capture a higher percentage of data variance after reducing the dimensionality. This is primarily because our method, unlike the PGA, does not require the submanifold to be a geodesic submanifold and to pass through the Fréchet mean of the data. Succinctly, the nested Grassmann structure allows us to fit the data to a larger class of submanifolds than PGA. We also proposed a supervised dimensionality reduction technique which simultaneously differentiates data classes while reducing dimensionality. Efficacy of our method was demonstrated on the OASIS Corpus Callosi data for dimensionality reduction and classification. We showed that our method outperforms the widely used PGA and the recently proposed PNSS by a large margin.